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/Abstract 

This report summarizes the attempts made to apply semidirect 
methods to the calculation of three-dimensional viscous flows over 
suction holes in laminar flow control surfaces. The attempts were 
all unsuccessful, due to either (1) lack of resolution capability, 
(2) lack of computer efficiency, or (3) instability. 
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INTRODUCTION 


Laminar flow control technology has the potential for signi- 
ficantly increasing the range and reducing the fuel consumption of 
commercial aircraft. With laminar flow control, the viscous drag 
of large surfaces such as wings is reduced by stabilizing the 
boundary layer, preventing or greatly delaying the transition from 

O' 

laminar to turbulent flow. 

On the conventional laminar flow airfoil, boundary layer sta- 
bility is enhanced by overall design of the surface so as to avoid 
pressure peaks and their inevitable concomitant, adverse pressure 
gradients. In the unconventional LFC technology, further stabili- 
zation is achieved by using a suction surface on the wing or fuse- 
lage surface to withdraw part of the boundary layer. The suction 
may be accomplished through long spanwise slots, a porous surface, 
or discrete cylindrical holes. 

The first-order prediction of the stabilization resulting 
from the surface suction is well understood by the classical 
boundary-layer stability theory based on the parallel-flow Orr- 
Sommerfeld equation, associated with C. :C. Lin, Tollmien, 
Schlichting, and others (e.g., Refs. 1, 2). Second-order effects, 
which could be significant for slots, may require the inclusion 
of x-dependency , as shown by Saric and Nayfeh (Ref. 3) . The dis- 
crete hole surface can be manufactured by electron-beam technology 
and has structural advantages over slots and economic advantages 
oyer^ porous surfaces manufactured by sintering techniques. 



However, the complicated three-dimensional flow over the cylindri- 
cal suction hole could possibly develop a secondary instability 
which would not be revealed by the boundary layer stability analy- 
sis of the unperturbed boundary layer flow. 

In order to study this secondary instability theoretically, 
an accurate representation of the three-dimensional flow in the 
neighborhood of the suction hole is required.- Particularly needed 
are accurate values of first and second partial derivatives in the 
y-direction normal to the surface and in the spanwise z-direction. 
Parametric studies require variation of the inflow boundary-layer 
thickness and profile, pressure gradients, hole diameter and span- 
wise spacing, skin thickness, and suction rates, as well as purely 
numerical parameters such as mesh-size. 

The objective of the present work was to develop a computer 
code to rapidly and accurately calculate the three-dimensional 
viscous flow over the suction holes, using the complete compres- 
sible Navier-Stokes equations. The intended technique was to be 
an application of the semidirect methods previously developed by 
the Principal Investigator. 


PREVIOUSLY USED METHODS 


The classical method of analysis for viscous flows involves 

the calculation of an inviscid flow followed by a boundary layer 
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calculation. This technique has been developed to a high level of 



sophistication in two dimensions and in three dimensions (e.g.. 
Refs. 5, 6) . With the inclusion of viscous iteration, it can even 
be used to calculate mildly separating flows in two dimensions 
(Ref. 7) . The boundary-layer method appears to be the method of 
choice for the suction slot problem. However, the discrete suc- 
tion hole introduces added ellipticity in the spanwise direction, 
and the boundary layer equations become less dependable and more 
difficult to calculate. 

The next alternative is to use the parabolic marching equa- 
tions... These equations are intermediate between the classical 
boundary layer equations and the full Navier-Stokes equations. 

This approach has been successfully applied to simple flows in 
constant-area ducts (Refs. 8, 9) and is quite economical compared 
to time-dependent Navier-Stokes solutions. However, there remains 
some questions as to the validity of the procedure in regard to 
the treatment of longitudinal pressure gradient terms. The prac- 
tical accuracy is suspect in the present case wherein those terms 
may vary significantly in the immediate neighborhood of the hole, 
which region is the most likely candidate for secondary instability 
(Ref. 4) . 

A third alternative is to use the full Navier-Stokes equations 
with no approximations, solved with a conventional time-dependent 
or time-like iteration scheme (e.g.. Ref. 10), obtaining the de- 
sired steady-state solution as the asymptotic limit of the time- 
dependent problem. However, three-dimensional time-dependent 



calculations of viscous flows are notoriously expensive, especially 
for the high resolution in y required for the present problem. 


SEMIDIRECT METHODS 

The present work involved "semidirect methods" for solving 
the full Navier-Stokes equations for the steady-state solution, 
using iterative procedures which are not time-dependent or even 
time-like. The principal investigator had previously developed 
several such methods with considerable success. (See Refs. 11-15. 
See also Ref. 16, and similar methods for the transonic inviscid 
equations in Refs. 17 and 18.) 

The basic concept in semidirect methods is to use the recently 
developed fast (direct, or non-iterative) linear equation solvers 
to solve linearized steady-state equations, which are then iterated 
to remove the non-linearities. This iteration is not time-like, 
as shown in Refs. 11-13. Convergence can be obtained in as few as 
6 to 8 iterations for some problems, and more typically in 10 to 
15 iterations for problems with flow-through computational boun- 
daries like the ones of interest here. Each of these iterations 
requires computer time comparable to that for a single time-step 
in an efficient time-dependent method, i.e. one which uses a direct 
elliptic solver for the Poisson equation. If the more common 
iterative methods (ADI, SOR, etc.) are used for the Poisson equa- 
tion in a time-dependent method, then the semidirect methods can 



sometimes attain a complete steady-state solution with computer 
time comparable to that of a single time step. 

The basic concept of the semidirect iteration scheme can be 
described by reference to a two-dimensional laminar flow problem 
using the vorticity transport equation, 

C t = -Re V*V£ + V 2 c 

where £ is vorticity. Re is the Reynolds number, and V is the vec- 
tor velocity. A steady state is assumed, so the time derivative 
is set =0. In the Split NOS method proposed here, the velo- 
city V at the k-th iteration is split into an initial guess V 
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and a perturbation V ' , not necessarily small. That is, 

v = V + V' 

Then the steady-state form of the vorticity transport equation is 
solved for the k-th iteration using 0(/Ax ) centered differences, 
as 

L(c k ) = Re 7 ,» (V'£) k-1 
where the linear operator L is defined by 

L (?) = V 2 C - Re V* (V°C) 

For the present 3D problem, the treatment of the spanwise convec- 
tion terms prevented a successful calculation, as will be described 


later. 



PARTICULARS OF THE 
LFC PROBLEM FORMULATION 

For the presently considered LFC problems, the free-stream 
pressure is prescribed, so that the primitive (velocity-pressure) 
system of variables is recommended. The principal investigator 
had also applied the Split NOS method to primitive variables in 
incompressible flow with success in 2D. The additional equation 
for energy required for the compressible flow solution does not 
present any special difficulties, since the boundary conditions 
(on heat flux or temperature) are well specified. (Free-convection 
problems have recently been calculated by the P.I. using semidirect 
methods.) The compressibility correction terms were to be lagged 
in the iterations, as are the velocity perturbation terms above. 

The work of Ref. 17 on inviscid flow indicates that compressibility 
effects in subsonic and slightly supersonic flows are compatible 
with a semidirect iteration scheme. 

Because the subject problem involved flow-through computation- 
al boundaries, rapid convergence requires (Refs. 11-15) that the 
linear solver used can accomodate first-derivative (convection) 
terms at least in the primary flow direction. Also, accurate re- 
solution of derivatives in y requires a fine mesh and possibly 
coordinate stretching in y. These requirements dictated that the 
commonly available fast Poisson solvers are not applicable to the 
proposed solution method. 

The EVP method previously developed by the Principal Investi- 
gator (Refs. 19, 20) and recently extended to three dimensions 



(Ref. 20-23) is applicable to this problem. Three different 3D 
linear solvers were utilized in the present investigation; simple 
3D marching, 3D marching with a banding approximation, and the 3D 
EVP/FFT method. It was intended that the solution to the title 
problem would be obtained with the full Navier-Stokes equations, 
rather than with some parabolized approximation; this will likely 
be significant near the suction hole. The constant-property 
approximation was used. This assumption is not necessary for the 
semidirect method, but the simplification helped the code debugging 
process . 

The equations used were the three momenta equations and the 
total energy equation, all in conservation form, and a modified 
Poisson equation for pressure. It is believed that the use of the 
pressure equation, rather than a steady-state continuity equation 
for density, would improve the accuracy for the slightly compres- 
sible flows of interest, would be compatible with the specification 
of the pressure in the external inviscid flow, and would circumvent 
the problem associated with the semidirect iteration of the con- 
tinuity equation as described in Ref. 12, page 194. The surface 
boundary conditions on temperature could be either specified 
temperature or heat transfer boundary conditions, including adia- 
batic walls. These and other boundary conditions are the same as 
would be used in a conventional time-dependent approach. In this 
regard, the conventional extrapolation conditions at outflow 
boundaries (e.g., Ref. 10) are the simplest to program and to 
adapt to a vectorized computer. 



The most recent work on semidirect methods for viscous flows 


(Ref. 15) indicated no numerical stability limitation on the 2D 
iteration scheme up to Re = 10 7 . The steady-state solutions 
obtained are the same as those obtained by the usual time-dependent 
methods using the same spatial differences, same variables, and 
same mesh. 

Iterative convergence rates would be improved (even in 2D) if 
all the terms possible were retained in the linear operators of 
the Split NOS method. However, it was decided to lag in the itera- 
tions the terms in the momentum equations which are the departure 
from incompressibility and constant viscosity, and those in the 
energy (temperature, not total energy) equation which arise from 
dissipation effects. This would result in some decreased itera- 
tive performance at high M, but would save considerably in storage. 
If all the possible terms were included in the linear operator, 5 
large matrices at each of the z-wave numbers would have to be 
stored and solved by LU decomposition, one for each elliptic equa- 
tion for u, v, w, T and P. With the ' scheme described above, the 
same matrix would serve for u, v, w and T, with a separate matrix 
for P. . Thus, only 2 matrices would serve for each of the five 
variables at each wave number. Also, the 3D EVP/FFT method itself 
has a savings of almost h because of the symmetry in the z-wave, 
number (Ref. 23) . Thus, the storage penalty for '.this method would 
only be about 20% of the basic storage required for u, v, w, T and 
P. This could be further reduced, by using Hockney's method for 
the pressure equation, to about 10%. For large problems to be run 



on the STAR 100, the savings in storage would translate into 
savings in time spent for page transfers, and would likely decrease 
actual computation time more than enough to compensate for the in- 
creased number of iterations. 

The choice for the mesh for the pressure equation was made. 
After considerable investigation, it was decided that a staggered 
grid for the pressure is very desirable, even though a non- 
staggered grid (in which all variables are defined at the same 
locations) is typically used in compressible flow solutions. These 
conventional solutions invariably utilize an explicit or an inher- 
ent artificial damping for the purpose of shock-capturing and 
stabilization. This artificial damping also depresses the parasi- 
tic mode that arises from the pressure solution. In the present 
case, the solutions would not contain artificial damping. How- 
ever, the MAC-type grid was also rejected because of its complex- 
ity in 3D, which is especially difficult for the interpolative 
coupling between the flow above the plate and the flow in the suc- 
tion hole. Also, the MAC grid would require different operators 
for u, v, w and T, with the resulting storage and time penalties 
mentioned above. Therefore, it was decided to use a grid in which 
only pressure is staggered, while u, v, w and T are all defined on 
the same grid. There is some penalty in operation count, but the 
interpolative coupling problem is tractable and the parasitic 
solution mode does not exist. 

It was intended to perform the calculations in several grids, 
which would be patched together either iteratively or directly. 



A cylindrical coordinate grid was to be used, inside the suction hole. 
Above the wing surface, a fine-grid and a coarse-grid solution would 
be obtained in rectangular coordinates . 

The use of the rectangular coordinate system above the wing 
surface would allow the streamlines to roughly follow coordinates, 
resulting in better accuracy and aiding (it was thought) conver- 
gence of the LAD-like iteration in z, as previously described. 

The several methods of solution attempted in this work will 
be described in detail in the following sections. The motivations 
and decisions about these methods are somewhat esoteric because 
of their intimate connection with the direct methods used for the 
linearized vequationst. The characteristics of these linear solvers 
are described in detail in Refs. 21-23. It should be realized 
that those methods are not easy to program in 3D, and that most 
of the P.I.'s time was spent in coding and debugging. 


3D MARCHING: 

FULL MATRIX AND BANDED APPROXIMATION 

The "simple” 3D marching method (Refs. 20, 23) was coded and 
tested. Previously, the 3D EVP method had been worked out theore- 
tically but had not been validated in an actual computation. 

The operation count for the 3D EVP method in its simplest 
form is very bad -with initialization for an MxMxM problem requiring 
0(M 6 ) operations. It would be necessary to use a banding approxi- , 
mation to the full influence coefficient matrix in 3D to make the 



method economical. This banding approximation had been proven in 
2D and worked out theoretically in 3D? it was now proven in 3D. 

The 2D experience carried over with no surprises to 3D, with only 
the usual programming difficulties. The result is the only 3D 
elliptic solver with optimal operation count, 0(M 3 ), for repeat 
solutions. Initialization requires 0(M 4 ) operations. 

The validating experiments with various bandwidths were first 
performed by simply zeroing the terms in the full influence coef- 
ficient matrix which were outside the band; i.e., the calculations 
were performed with a full matrix solution routine with null cal- 
culations. Later the code was actually converted to a banded 
matrix solver to realize the advantages in operation count and 
storage. The relative advantages of using a simple banded matrix 
solver vs. a block-banded solver were studied. The former appeared 
to be preferable. The available banded matrix solvers were studied 
and the IMSL programs were rejected. in favor of the LINPACK pro- 
grams written by Prof. Cleve Moler of the University of New Mexico 
Math Department. His programs give an estimate of the condition 
number of the matrix problem, and are well-designed to avoid 
unnecessary accumulation of round-off error, which is critical 
in the large 3D problems. 

After considerable difficulty with operating systems, the 
LINPACK system was set up on both the NOS system and the batch 
CDC 6600 system. The latter proved to be necessary because of 
the stringent storage limitation on the NOS system. 



A large parametric study was then performed of the accuracy 
of the solutions. Both the condition number of the matrices and 
the maximum error in the field were investigated as functions of 
cell aspect ratio, independent mesh dimentions I, J and K, band 
width, number of corrective iterations, and dimensionality. Also 
studied were the solution time for initialization, solution time 
for repeats, and storage. 

The adequacy of the banding approximation is much more sensi- 
tive to J, the problem size in the marching direction, then had 
been previously realized. Also, with KBAND = 1 (a minimal appro- 
ximation which includes only the nearest neighbors in the banding 
approximation) a 31x31x31 problem does not fit into core in a 
CDC 6600, since the storage penalty is about 3*KBAND*M 3 . A 
26x26x26 problem fits, but .the condition of the matrix is so poor, 
even at cell aspect ratios = 10, that accurate solutions are not 
possible. The best problem that fits is 21x21x21 with KBAND = 3, 
which gives excellent accuracy of the finite difference solution. 
However, it does not seem that this resolution would be adequate 
for the LFC problem. 

A scheme was then analyzed for effectively doubling the mesh 
size by solving on a submesh and using iterative corrections. The 
method appears to be useable, but doubles the mesh resolution only 
in one direction. It would also be very difficult to code. Fur- 
ther, a discussion on 2/16/78 with D. Bushnell indicated that an 
OCA 1 *) solution even in a 31x31x41 mesh, while adequate for initial 
studies, would hot be adequate for final calculations. The methods 



considered are mesh-size limited, i.e. , the submesh scheme cannot 


be applied recursively in 3D. It thus became evident that high 
mesh resolution with the semidirect methods in 3D was going to 
be very difficult. High resolution is not a difficulty in 2D, and 
the codes in 2D can be arranged to systematically refine the mesh, 
but the techniques are not simply extendable to 3D. Thus, the 3D 
solutions require different nonlinear iteration methods dictated 
by the linear solvers. . 


3D EVP/FFT 
LINEAR SOLVER 

The remaining nonlinear iterative methods which were attempted 
used the 3D EVP/FFT linear solver. The peculiar properties of this 
solver limit the success of the nonlinear iterations, so it is 
necessary to describe. its characteristics here. 

The 3D EVP/FFT method is a combination of the original two- 
dimensional EVP method (Ref. 19) and Hockney's method (Ref. 24). 

The dependent variable is first Fourier transformed in the z dir- 
ection using the Fast Fourier Transform or FFT. Then the EVP 
method is used to solve the two-dimensional problem in x and y 
for each Fourier component, followed by a reverse transformation. 

The operation count for this method for an IxJxK problem is 
0(J • J *K • InK) for repeats. Although the Fourier-transformed var- 
iable is complex, the required influence coefficient matrix Cl 
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(Ref. 23) is real. Further, it is symmetric in the wave number 
k , thus reducing the required storage and operation count by 
almost a factor of The three-dimensional equations require 

storage and initialization of only two Cl's, one for the three 
velocity components and one for the pressure. The error propoga- 
tion characteristics of this three dimensional problem were worked 
out in detail (Ref. 23) and found to be applicable to the coordin- 
ate and mesh systems of the present problem. This EVP/FFT direct 
linear solver retains the ability of the EVP method to handle 
variable-coefficient first-'.and cross-derivative terms in the x-y 
plane, although it is restricted to constant-coefficient second- 
derivatives in z, as are the more common methods. 


TIME-LIKE AND LOCALLY-ONE- 
-DIMENS IONAL ' METHODS 

In the interest of solving the LFC problem as efficiently as 
possible, several schemes were considered which fell somewhere 
between the semidirect methods and the more conventional time- 
like methods. 

An obvious approach is to use 2D methods in x-y, iteratively 
coupled to the terms in z, by way of a tridiagonal or hopscotch 
solution in z. However, there is no reason to expect such a 
planar ADI or planar hopscotch method to converge any faster than 
a 2D ADI solution, so this approach was discarded, even though 



the ease of coding was attractive. Another method which would be 
easy to code is a BIR (block-implicit relaxation) method extended 
to 3D. Here, the requirement for high resolution would again slow 
convergence, so this approach was discarded. Several variants of 
a locally-one-dimensional correction scheme were devised and tested. 
The first scheme devised was an iteration scheme which can be de- 
rived as a LOD method with an infinite time step. The first h~ 
step uses the 3D EVP/FFT linear solution. Then, instead of merely 
lagging all the cross flow advection terms, these are used in the 
next %-step in a ID correction which includes advection. Part of 
the z-diffusion term is needed in the second %-step in order for 
that term to be elliptic. The method was programmed and tested 
with a general weighting factor W in the usual manner, with 0 <_ 

W <_ 1, to split the z-diffusion terms between the first and second 
half-steps. The code also allowed for finite At, for different 
At's in the first and second %-steps, and for different numbers 
of cycles in the first and second %-steps. For example, the code 
accepts a single cycle at At^ for the first %-step, and n cycles 
at At 2 = At^/2 for the second %-step. The experimentation was 
performed on a linear constant coefficient 3D advection-dif fusion 
equation, with input parameters of directional Reynolds numbers, 

Rx , Ry , and Rz . 

The results were uniformly negative. Extensive parametric 
tests were run without discovering a stable combination, even for 
a li- Reynolds numbers = 0. This was surprising, and there is always 
the possibility of coding errors, but the store of debugging 



techniques and much of the computing budget was exhausted, and the 
work proceeded on the premise that these LOD methods are in fact 
unstable. 


SEMIDIRECT ITERATION WITH 
LAGGING SPANWISE - ADVECTION TERMS 

The final method which was coded and tested was a combination 
of the Split NOS and LAD methods. The linear solver used was the 
3D EVP/FFT method, which allows only 3 2 /3z terms in z, and limits 
the advection coefficients in the linear operator to u(x,y) and 
v(x,y) . In earlier published results on the LAD method in 2D i 
(which contains no advection terms in the linear operator) conver- 
gence was limited to low Re. However, this method was analyzed 
(under separate funding) and modified to achieve convergence at 
high Re, although at a loss in convergence speed. 

The 3D EVP/FFT routine was coded into a CALLable subroutine, 
and used as the basis of a fully 3D fluid dynamics solution proce- 
dure using the velocity and pressure variables. Extensive experi- 
mentation was carried out with this code and with a related 3D 
code, developed under separate contract, using the vorticity- 
velocity variables. With both codes, tests were conducted on the 
simple problem of perturbed 3D Couette flow. 

It became evident that the 3D EVP/FFT method would not be 
useable for 3D fluid dynamics solutions using semidirect methods. 



Instability resulted at all non-zero Reynolds numbers in both 
codes, using a variety of boundary conditions. When the pressure 
solution was removed from the nonlinear iteration procedure, sta- 
bility was slightly improved, but still not useable. With only 
the u-component of velocity active, stability was achieved only 
through Re = 0(10) . For Re = 15, convergence was obtained but 
very slowly and erratically. There always remains the possibility 
that the failure was due to coding errors, but the test for only 
the u-component is fairly straightforward; i.e., there were no 
changes in boundary conditions (except when periodic conditions 
were used in the third direction) and no interactions with other 
velocity components or pressure. When the equation was linearized, 
correct solutions resulted. Also, in 2D, the method is convergent 
even for high Re, for the nonlinear equations. It therefore 
appears that the instability is inherent in the 3D solution using 
the 3D EVP/FFT direct solver. It is conjectured at this time that 
the instability is related to the sensitivity of the FFT to small 
disturbances . 


OTHER LFC APPLICATIONS 


The semidirect solution methods still have these possible 
applications to the LFC suction boundary-layer control problem. 

(1) 2D solutions can be obtained very rapidly, and with high 
resolution. (2) The linear 3D EVP/FFT method can be used for the 



direct solution of the pressure field for use with a time-dependent 
or time-like solution procedure in 3D. 
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